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ABSTRACT 


The role of non-destructive testing { NDT ) in high technology 
industry is crucial and inescapable. The quality consciousness of the 
structural products in both the manufacturing stage and during the 
operating lifetime are the main reasons for the industry to be aware 
of this (NDT) tool. Two of the key areas where NDT is pivotal are 
aeronautics and nuclear research. Ultrasonic means have found to be 
very effective in the inspection of composite materials and austenitic 
steels, both are highly anisotropic in nature. Scattering know^ledge of 
ultrasonic waves by various flaws provides all the information 
regarding the inverse characterization of the defects which has 
not been thoroughly solved yet. 

Computer simulation is a powerful aid to the understanding of 
physical processes. In the present work it is shown how Finite Element 
Analysis can be used in the study of ultrasonic wave propagation and 
scattering in general anisotropic media. Pulsed compression and 
Rayleigh wave interactions with surface cracks and body defects in 
isotropic media has been modelled. Also studied is the scattering by a 
delamination and a hole in a orthotropic composite laminate. These 
simulations are based on the numerical approximation of the hyperbolic 
equations of infinitesimal elasticity with simplest boundary 
conditions. The two-dimensional FEM model takes account of short 
duration pulses, finite-width beams and of different polarisations. 

The physics of scattering is clear from the results, which are 
shown in the form of snapshots to visualise the mode-conversion, 
reflection and diffraction of scattered ultrasonic waves. The 
qualitative results are in agreement with the photo elastic fringe 
patterns and experimental observations. A marked difference in the 
lobe structure of scattering cross sections of composite laminates 
explains the beam skewing effect in anisotropic media. - 
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CHAPTER I 


INTRODUCTION 


1.1 GENERAL INTRODUCTION 


The elastic wave propagation and 
investigation for several years, because 
associated with the seismology, nuclear 
response to high rates of loading, and 
The effects of very short term impact 


scattering have been under 
of the technological needs 
test monitoring, structural 
non destructive evaluation, 
loads on structures can be 


understood fully if analyzed in terms of propagation of stress waves. 
Crack propagation or the interaction of dynamic stress fields with 
existing cracks, voids or inclusions in a material can be predicted if 
the dynamic effects are taken into account. Further the knowledge of 
elastic wave propagation is essential in a host of other fields such as 
ultrasonics, geophysics, electronics etc. 


The reliability and suitability/integrity of a structural component 
for its further service depends upon the current material state, failure 
modes and service conditions. Cracks are one feature which limit the 
strength of structures under stress and often causes its catastrophic 
failure. But voids, inclusions, interfaces, distribution of cracks and 
in general damaged regions of material will have dangerous effects on 
the performance of the components. Critical crack sizes vary from 
several tens of millimeters in thick steel structures dowm to fractions 
of a millimeter in turbine blades and other aero components, while in 
the case of ceramic materials these are of the order of few micrometers. 
Thus it is highly desirable to assess the material state. The role of 
quantitative non destructive measurements lies in identifying and sizing 
microscopic flaws that would ultimately lead to failure. 

Because of their ability to penetrate into the interior of the 
metal, composite material or ceramic parts the elastic waves are one of 
the most important forms of probing energy. Ultrasonic waves are such 
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types ol waves which exhibit rich propagation phenomena and can be 
utilized in extracting a wide range of information about the structural 
integrity ot a component. The scattering, reflection and mode conversion 
of are the different phenomenon of pulsed elastic waves that provide 
the basis for ultrasonic NDT. 


In genei'aJ , ultrasonic flaw detection and characterization 
techniques can be carried in two approaches. The imaging approach seeks 
to process the scattered field in such a manner that a visual outline of 
the object is produced on a display. The inverse-scattering approach 
attempts to infer geometrical characteristics of a flaw from either the 
angular dependence of its far-field scattering amplitude at a fixed 
frequency, or from the frequency dependence of its far-field scattering 
amplitude at fixed angles. The solution to the direct scattering 
problem, that is, the computation of the field generated when an ultra 
sound wave is scattered by a known flaw, is a necessary step to 
understand the solution of an inverse problem. 

In most of the ultrasonic non-destructive applications the material 
is taken to be homogeneous except for any defect. The wave propagation 
is then rectilinear with simple laws of reflection at plane boundaries 
except for the complex scattering phenomenon at the inhomogeneity. If 
the grains are larger than an ultrasonic wave length and there is a 
preferred texture then the material becomes anisotropic. Examples of 
materials w^here anisotropy is apparent are composite materials such as 
fiber-reinforced plastics, austenitic steels, crystalline semiconductor 
materials and rocks. 


Composite materials are widely used in structural applications 
where high strength- to-weight and stif f ness-to-weight ratios are 
required. Air crafts are typical weight sensitive structures in which 
composite materials are cost effective. Composite materials can be 
tailored to efficiently meet the design requirements of stiffness, 
strength and many other parameters in various directions. 
Unfortunately, the manufacturing process of composites may result in the 
presence of unwanted artifacts defects such as voids, resin rich areas, 
distorted fibers and inclusions. Damage and general material degradation 
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can also occur during the in-service operation of composite components. 
Typical causes ol such damage are continuous cyciic/f atigue loading | 
rapid changes in local temperature and moisture, and impact loading 11]. 
So, tests are needed in both the manufacturing phase of composite 
structures in order to assure quality as well as during the operating 
lifetime ot the components in order to predict and hence prevent 
failure . 

Propagation of waves through homogeneous anisotropic media is more 
complicated than propagation in homogeneous isotropic media. The 
direction of energy flow is not necessarily at right angles to the wave 
front. This effect, known as beam skewing^ has been interpreted as 
causing much of the difficulty in the inspection of composite 
structures. A large number of wave propagation studies in a laminated 
and fibrous composite medium are limited by the assumption that the 
structure is periodic or the distribution is detei-minist ic [2J. 

The motivation of the present work is to study the w^ave propagation 
and scattering of ultrasonic waves in materials which are both 
anisotropic and inhomogeneous but v/ith a scale of inhomogeneity varying 
from a fraction to a few times ' a ultrasonic wavelength. The defects in 
composite materials are just of this nature. 

1.2 ULTRASONIC INSPECTION 

In ultrasonic inspection short bursts of a few cycles of elastic 
wave energy, typically in the frequency range of 1-10 M Hz, are injected 
into the part being inspected. The variation in the speed of wave 
propagation and energy loss by interaction with material microstructure 
are the w^ave parameters obtained for the characterization of material 
properties. A typical Ultrasonic NDT system is illustrated in Fig 1.1. 
Flaws are indicated by the reflection of the energy (echoes) back to the 
same transducer (pulse-echo) as shown or to a second transducer 
( pitch-catch/ through- transmission ) . In state-of-the-art sy^stems, the 
signal processing takes the form of amplification followed by 
rectification L3J. Flaws are detected when an echo appears that exceeds 
a preselected threshold. The size is then estimated based on either the 
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peak value ol this echo or the variation of echo height as the probe is 
mo\ ed . 


Uni ortunately , quantitative characterization of a flaw i.e., the 
measurement ol size, location, orientation and roughness, based on 
amplitude measurements may be v’ery inaccurate unless the measurement is 
careiully calibrated. The reason is that part of the incident energy 
that insoniiies the ilaw is scattered rather than specularly reflected. 
Hence, the calculation of the total field which includes scattering 
eiiects is rather complicated. It was also an established fact that the 
attenuation coetiicient of an ultrasonic wave highly depends upon 
scattering rather than absorption (due to internal friction and 
dislocation movement) at frequencies higher than 1 M Hz. 

Three elements play a major role in the physical solutions of 
detecting and sizing flaws : transducer improv'ements , high-speed digital 
processing and elastic wave scattering theories. In practice, 
transducers do not follow theoretical expectations 14 J. Hence, their 
individual properties (radiation pattern, impulse response etc.,) must 
be measured and taken into account when interpreting experimental data. 

The essential feature of improved interpretation of ultrasonic 
signals is a fundamental understanding of elastic wave - flaw 
interaction. The limitations imposed by the analytical methods in these 
studies have been overcome in the recent years by numerical methods, 
which require high speed digital computing facilities. The whole range 
of modelling techniques are becoming increasingly attractive and 
practical with the dev’^elopments in and the availability of super 
computers. Modelling studies can play an important role in Ultrasonic 
NDT by providing : 
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Fig 1.2 Model of general NDT 
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TABLE 1 


Range ot wavelengths and order of magnitude for feature size in studies 
of wave propagation . and scattering in solids 


Subject 

Wave Length 
(X) 

Frequency 

(f ) 

Feature 
Size (d) 

Order of 
Magnitude { M ) 

Geophysics 

Seismology 

4 0 km 

'■ 0.05 Hz 

1000 km 
10-100 m 

10^ m 

10" m 

Civil Eng. 
Mech. Eng. 

10 m 

100s Hz 

1-10 m 

10 m 

NDT 

0 . 1-100 mm 

0.025-15 MHz 

1-10 mm 

10'- m 

Electronics 
( SAW ) 

0.05 mm 

>40 MHz 

0.01 mm 

10‘- m 


ULTRASONIC NDT/NDE 

Defect characterization 
Material property studies 
Acoustic emission 
Acoustic microscopy 



Site investigation 

Vibration studies 

Earth quake damage prediction 


Fundamental elastic wave 
propagation and scattering 
Fracture mechanics & 
Dynamic crack propagation 


Fig. 1.3 Fields of study which include the investigation of 
elastic wave propagation phenomena. 




1.3 LITERATURE REVIEW 


Understanding the propagation and scattering of waves through an 
elastic medium is of potential benefit in several fields. Fig 1.2 shows 
the fields of study which includes the investigation of elastic wave 
propagation phenomenon. An excellent research literature on the state of 
the art of stress wave propagation and scattering was reviewed by Pao 
[5J. The role of elastodynamic scattering problems in quantitative 
non-destructive evaluation has been reviewed by Gubernatis [6] and Bond 
et al 17 j. Some practical applications have been discussed by Coffey and 
Chappman [8 J . 

Most of the classical solutions of scattering by a flaw are 
applicable either at very high frequency range, 4a » 10, or at very low 
frequency, 4a « 1, where ‘4’ is the wave number and 'a’ is the flaw 
major dimension. In the mid-frequency range, where the wave length of 
the elastic wave is approximately equal to the flaw dimension, the 
solution of elastic scattering problem become very complex because of 
mode-conversion phenomenon involved. In addition, the analytical 
solutions are restricted to infinite or semi-infinite medium, 
homogeneous and isotropic solids and idealized defect shapes. Results 
have been obtained for wave interactions with a penny-shaped crack, 
infinite strip, cylinder, sphere and wedge [43,44] . 

The following are the commonly used techniques to solve various 
elastodynamic scattering problems. 

1.3.1 Separation of Variables 

The scattering of elastic waves by a circular cylinder or a sphere 
can be analyzed by expanding the scattered field in a series of wave 
functions [9J. For these shapes of flaws the wave equations are easily 
separable in an orthogonal curvilinear co-ordinate system to allow for 
the matching of boundary conditions. The solutions thus obtained have 
wide applications in analyzing dynamic stress concentration and 
interaction of structures with ground waves. 
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1.3.2 Volume Integrals 


The waves scattered by an inhomogeneity are expressed by a volume 
integral of equivalent body force multiplied by the Green’s function for 
the suri'ounding medium. The volume integral which contains the unknown 
total wave field can be solved by iteration, and the solution is called 
Born approximation, if the incident wave field is substituted for the 
unknown total wave as the first trial [10]. 

1.3.3 Surface Integrals 

The surface integral formulation is useful for treating crack-like 
geometries. The governing differential equation for elastic waves is 
converted to an integral over the surface of the scatterer. The surface 
integral equation is solved numerically by using Hilbert-Schmidt method 

[11] . It can also be reduced to a system of linear algebraic equations 
for the scattered amplitude which are then solved numerically by matrix 
methods. Two such methods are transition matrix (T-matrix) theory 

[12] and the method of optimal truncation (Moot) [13]. 

1.3.4 Diffraction Methods 

At high frequencies, the scattering of elastic waves are analyzed 
on the basis of elastodynamic ray theory [14]. In Keller’s [15] 
Geometric Theory of Diffraction (GTD), the diffracted field is related 
to the incident field through diffraction coefficients. These are 
determined through the solution of an appropriate canonical problem. 
However, the lack of canonical solutions for elastic waves imposes a 
serious limitation of GTD method. A fully three-dimensional ray tracing 
computer code for anisotropic and inhomogeneous media was developed by 
Oglivy [16] . The model takes into account correctly the amplitude and 
phase changes at boundaries between media and can therefore, be used to 
make quantitative as well as qualitative predictions. Ray tracing work 
has shown that angled compression and shear horizontal waves maintain 
their focus well, whereas shear vertical polarized waves are always 
greatly defocussed. 
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1.3.5 Numerical Methods 


These are the computer-modelling techniques which are based on the 
numerical solution of the partial differential equations describing the 
scattering system. In fact, these are the only methods available in the 
complex mid frequency region, where the scatterer dimensions are of the 
same order as the ultrasonic wavelength. A surv'ey of numerical methods 
used for solving the elastodynamics was discussed by Bond [17]; in 
particular to seismology was given by Chin et ai ; IIBJ. The following 
three families of techniques are commonly employed; 

1. Boundary Element Methods 

The Boundary Element Method involves the solution of an integral 
equation, which is obtained by dividing the surface of scatterer into 
small elements each of which then has an unknown displacement L19J. 
Schafbush et al [20J have applied three dimensional boundary element 
method to elastodynamic problems related to ultrasonic numerical 
differential equations. Hail and Robertson [21] have studied the 
scattering of a plane wave by a circular cylinder using 2D boundary 
element method and compared the results obtained for various elements 
with the analytical solution. 

2. Finite Difference Methods 

In numerical modelling applied to elastic waves much effort has 
been concentrated on the finite difference method. Finite difference 
models give a complete description of the time development of the 
wave-field of interest. The methodology involves the replacement of the 
governing partial differential equations by difference equations, and 
then solving them numerically. Alterman and Karal [22]; and Bore [23] 
have employed this method for seismology applications. 

Developments in the finite difference modelling of ultrasonic wave 
propagation have been reported by C.J. Bond and his co-workers [17]. 
They have modelled a range of compressional , shear and Rayleigh wave 
systems on isotropic, homogeneous media, involving surface slots and 
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cracks. The results are shown in the form of a display in which vector 
displacements are scaled to appear as lines in the direction of 
displacement, and follow the progress of an elastic wave through the 
medium. These results have been found to agree w’ith the photo-elastic 
visu lisations. In addition, immersion inspection systems have been 
modelled where tliere is propagation in both an elastic solid and a 
viscous fluid 124J. 

To reduce the problems associated with the use of limited extent 
spatial grids absorbing boundaries have been investigated and 
incorporated into the models 125J. 

Bond [26] has emphasized on the "mode conversion" and " total-wave" 
concepts in particular to ultrasonic NDT and has giv^en details of 
experimental results of crack characterisation by making use of the 
mode-conversion phenomenon. 

Temple [27] has used a three dimensional finite difference model to 
study the propagation and scattering of elastic waves of different 
polarizations from a cylindrical hole in a transversely austenitic steel 
w’elds . Results in the form of A-scan display and scattering 

cross-section are given and compared with the experimental results of 
White [ 28 j . 

Harumi and his co-workers [29] have produced various ‘lumped mass’ 
schemes w'hich appear to overcome many of the limitations on the 
scatterer shape imposed by the conventional cartesian co-ordinates. 
Oglivy and Temple [30] have applied both ray-tracing and finite 
difference methods for predicting the behaviour of ultrasound in typical 
manual-metal-arc-w’eldments . They have demonstrated that the two 

approaches are consistent and complementary. 

3. Finite Element Methods 

The superiority of finite element method lies in its ability to 
handle complex domains, defect shapes, inhomogeneities, material 
non-linearities, applying the boundary conditions etc. 
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Smith [31J has developed a finite element code to study the body 
wave propagation problems in seismology* He modelled a deep earthquake 
zone realistically in two dimensions and earthquakes simulated at depth. 

Lord and his co-workers have produced a scheme to solve the 
elastodynamic problem using a finite element formulation for spatial 
discretization and finite difference formulation for time marching. 
Ludwig and Lord [32] has reported the details of 2-D finite element code 
employing Newmark’s algorithm to evaluate time integration. Results 
were given in the form of simulated A-scan and visualisation of 
wave-fields for a longitudinal pulse at different instants. Transmitter 
and receiver models have also been incorporated to take care of 
transducer effects [33]. The results were comparable with the practical 
ultrasonic NDT measurements. You and Lord [34] have reported 
improvements in the scheme from computational point of view, by 
employing a central-difference algorithm for time integration and a 
"lumped mass matrix" to avoid matrix inversion. They have also studied 
wave interactions with a rectangular slot and surface breaking cracks in 
an isotropic medium. 

A detailed analysis of various computational aspects of different 
algorithms used in the finite element method for the wave equation, from 
the accuracy, speed and storage points of view, has been reported by 
Seron et al [35]. It is concluded that, although central-difference 
scheme is the best in terms of computational speed and cost, Nev/mark’s 
constant average acceleration (CAA) scheme along with quadrilateral mesh 
gives the best accuracy. 

Ludwig et al [36] have reported a 2-D finite element modelling of 
wave propagation in an orthotropic material. Results were discussed with 
respect to the focussing, defocussing and skewing of ultrasonic beam. 

You and Lord [37] have demonstrated a finite element test bed for 
diffraction Tomography as an inverse-algorithm and used the normal 
component of the displacement data on the measure line to reconstruct an 
object function. 
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Recently Jaieel et al 13BJ have developed a Finite Element code to 
study the pulsed compression wave and Rayleigh wave interaction with 
cracks in isotropic and Graphite/epoxy composite laminate . It was found 
that for quasi-longitudinal wave the displacement along the fiber is 
prominent, whereas for quasi-transverse wave the displacement 
perpendicular to the fiber is prominent. 

An important parameter in elastic wave scattering is the total 
power scattered by the flaw. It is customary to relate the scattered 
power to the power per unit area normal to the in coming wave. This 
ratio has the dimension of area and is known as scattering 
cross-section, Barrat and Collins [39] have considered the total field 
as a superposition of a plane wave field and scattered field. They have 
show'n that cross section can be calculated from the far tield amplitude 
of the scattered wav'e in the direction of propagation of the incident 
wave. Krenk and Schmidt [40] have reported an extensive results of 
elastic wave scattered by a circular crack in the form of maps of the 
scattered far field as scattering cross sections. White[29] has 
calculated the cross sections for the case of scattering by a hollow 
circular cylinder, fluid filled cylinder and the results are compared 
with his experiments. 

Su et al. [41] have extended the unimomemt method to elatic wave 
scattering problems. This method uses the finite element approximations 
for the field in the vicinity of the scatterer and the scattered field 
outside this region by wave function expansions that satisfy appropriate 
radiation conditions at infinity. Goetschel et al [42] have employed 
exterior wave functions to the wave field with finite element to the 
near field to study the cross sections for a variety of axisymmectric 
inclusions . 

1.4 PRESENT WORK AND THESIS LAYOUT 

In the present work a two dimensional plane strain finite element 
model is developed to study the ultrasonic wave propagation and 
scattering in linear, elastic medium. 
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The model is capable of predicting the propagation of short pulses 
of ultrasonic waves in a finite width beam through a anisotropic medium 
and of predicting the amplitudes of waves scattered from wedge crack, 
side drilled holes type defects. 

Pulsed compression and shear wave interaction with a wedge crack 
and cylindrical hole in an isotropic medium were studied. The scattered 
waves are clearly showing the mode conversion phenomenon. Elastic wave 
propagation and scattering in Graphite/Epoxy composite has been 
investigated. The material is modelled as homogeneous and orthotropic. 
Differential scattering cross section were drawn for all the above 
mentioned cases. 

Chapter II of the thesis gives the basics of ultrasonic wave 
propagation and scattering. Wave propagation iii elastic medium, plane 
w’aves in isotropic and transversely isotropic medium were discussed. 
Formulation of boundary value problem is given. Ultrasonic wave 
scattering by a single flaw is discussed and scattering cross section 
details were also described. 

In Chapter III the finite element formulation of the hyperbolic 
wave equation is presented. Transducer considerations and artificial 
boundary problem are mentioned. Methodology employed to calculate the 
scattering cross section and computer implementation are also discussed. 

Chapter IV gives the results and discussions of the simulations 
performed using the FEM code. Conclusions and scope for future work are 
given in the last chapter V. 
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CHAPTER II 


BASICS OF PROPAGATION AND SCATTERING 
OF ULTRASONIC WAVES 


When a solid is subjected to dynamic loads, the stresses become 
functions of both space and time, and the elastodynamics is that branch 
of solid mechanics, which determine these stresses. The mathematical 
treatment of elastodynamic problems is distinguished from that of the 
elastostatic case by the presence of the the acceleration term 
Cauchy’s equation of motion [431. 


2.1 FIELD EQUATIONS OF LINEAR ELASTODYNAMICS 

Under the assumption of linearity, for a body of volume V enclosed 

by surface S = S + S , where on S tractions are prescribed and on 

g u (7 

displacements are prescribed, the equation of motion 144] may be 
written as 

( 2 . 1 ) 


a . . .+pb. =:pu. 

ij , J 1 1 » 


tt 


W^here a . is the stress tensor, u. is the displacement vector, b^ 
1 J ^ 

is the body force per unit volume and p is the mass density of the 
material . 

The linear constitute equations for a hyperelastic solid in the 
absence of thermal and other non-mechanical effects is of the form 


^ijki ^kl ^ijki 


( 2 . 2 ) 


Where C 
properties 


i iki 


is the elastic stiffness tensor with the symmetry 


c = c ~ c = c 

ijkl ~ jikl 


(2.3 ) 
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upon substituting (Z) in (1) 


a 


dx 

j 


C 

i jk i 




Sx 


+ p . 




at 


(2.4) 


This is the Navier-Cauchy equation of motion for an anisotropic and 
inhomogeneous elastic solid. It will not be possible to seek an 
analytical solution to these three coupled equations. Solutions in 
directions of high symmetry can be found for homogeneous media such as 
single crystals 145J, but even for arbitrary directions in homogeneous 
anisotropic media the solutions are best found numerically. 


For an isotropic solid, the elastic stiffness tensor has the 
special form 


i i k 1 


X 6 . 6, . + p f ^ 1 ^ ivA 

k 1 \ 1 k 1 1 i 1 J k j 

-j— - j 


(2.5 ) 


Where X, fj are the Lame’s constants and the stress strain relation 
becomes 

(2.6) 


a = X e, ,6. . + 2 fj €. , 

i 1 k k 1 1 1 j 


The strain displacement relations are 

e =4- r u. . . u. . 1 

i J ^ (^ 1 > J J > ^ J 


(2.7 


Substitution from (2.2) and (2.7) into (2.1) leads to the governing 
equation of motion in terms of displacements as 


X + /J u. ,. + fJ u + P 

[ J J > Ji 1 » JJ ‘ 


= P 


a-u. 


a t‘ 


( 2 . 8 ) 
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1 he vector equivalent of this expression is 

I X + 2/j I VV.lt - p VxVxTJ + p = P ' - (2. 9) 

" ^ at- 

The equation of motion {Z*9) represents a system of three coupled 
partial diiferential equations for the displacement components; they 
must be solved subjected to boundary conditions, 

. n. = t, on (Traction B.C.) 

- u] on S (Displacement B.C.) (2.10) 

- 1 u 

and initial conditions, 

u (x,U) = u” (x,0) 

i i 

u(x,0)=u°(x,0) (2.11) 

I i 

where t. , u’ , u° and u° are prescribed quantities. 

! 1 i i 

By defining the dilation and rotation vectors as 

A = V.u = €. . 

1 1 

1 _ 

0 ) = — r— Vxu 

The equation of motion (2.9) may be expressed as 

[ X + 2/J ]vA - 2p Axt!^ + p = p ^ (2.14) 

t J at 2 

The above form of representation explicitly displays the dilation 
and rotation. Another advantage is that the result is valid in any 
curvilinear coordinate system. 


( 2 . 12 ) 

(2.13) 


17 



2 . 2 FORMULATION OF BOUNDARY-VALUE PROBLEM 


The direct solution of the governing elastodynamic equations (2.y) 
IS quite dilficult due to the strong coupling between the three 
components ol the displacement vector u. A convenient way for 
decoupling them is to introduce scalar and vector potentials , by using 
Helmholtz's theorem. The displacement and body force vector can be 
wi'itten as 

u = grad ip + curl ip (2.15) 

where V. tp = 0 . In the above equation <p and ip are scalar and vector 
valued functions of position and time. 

F = grad 0 + curl x (2.16) 

The completeness of this representation has been discussed by Achenbach 
144 J. Substitution of u and F in (2.9) and after simplifying the terms 
leads to the following equation 

p 

grad I" c^ 'P- + 0 ] + curl f c^ . +,(1=0 (2.17) 

'' ■ at- V - at- > 

where c and c are longitudinal and shear wave velocities. 

1 2 

Clearly, if 


• • 

(P - c- v~<p = 0 (2.18) 

1 

and 

V - c^ y-jj) = X (2.19) 

then u satisfies the equation of motion. 

Hence <p and ip satisfy wave equations with two distinct wave speeds 
c^ and c^ respectively. Since c^ > c^ , these are referred to as 
P-( primary) waves and S- { secondary ) waves in geophysics literature. In 
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mechanics literature these are known as pressure or diiationai or 
compressionai waves and shear or equi voluminal waves, respectively. It 
is worth mentioning that the mechanical disturbances associated with 
these two tields 0 and ip propagate unchanged only if there are no 
boundaries or scattering surfaces present i.e., anisotropic and 
homogeneous medium. Boundary conditions at such interfaces will, in 
general, couple the two fields, giving rise to the phenomenon of mode 
conversion. In other words, either P-waves or S-waves will generate both 
P-waves and S“waves on striking a boundary. This is a feature which does 
not have a counterpart in acoustics and electromagnetics. 


2.2.1 Two-Dimensional Problems 


If the field variables u,v,w are independent of one of the 
coordinates, then the equations of elastodynamics can be split into two 
independent systems of differential equations and boundary conditions; 
one system describes the in-plane deformation and the other, the 
antiplane deformation of the solid. Assuming that the displacement and 
stress components are independent of the coordinate z, then the 
components of the Navier^s equation for an isotropic solid becomes, 




a A 


1 >. + H I 


I" X + H ] 


^ ax 

a A 


+ /J v~u = p 


+ ^ v-^v = p 


^ ay 


a-u 


at' 


a~v 
at 2 




= P 


a-w 

at- 


where A = u + v 
» X > y 


(2. 20. a) 


(2.20.b} 

(2.20.C) 
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Fig.2..1 Coordinate systems: (a) Cartesian coordinates; 

(b) two-dimensional XY (P-SV) coordinates; 

(c) two-dimensional XZ (SH) coordinates. 


2 


The in~piane displacement components u and v satisfy the first two 
partial ditferential equations and the in-plane stress components. Thus, 
it the boundary conditions of the problem are given in terms of these 
variables, then the displacement component will play no role in the 
solution ot the problem. In seismological problems, the xz-plane is 
otten taken as the horizontal plane on the surface of the earth. Since 
the displacement vector { u , v , 0 ) lies in the vertical plane through the 
direction of propagation, the S waves represented by the vector 
potential ip are called vertically polarized shear waves, or, simply SV 
waves . 


The antiplane displacement w satisfies the third partial 
differential equation ( 20 .c) and the antiplane stress components a and 

X 7 

^^ 7 ' Since the displacement component w is horizontal and is 
perpendicular to the direction of propagation, the waves represented by 
w are called horizontally polarized shear waves, or SH waves. Fig. 2.1 
shows various combinations of displacements in a cartesian coordinate 
system, which are generally employed. 

2.3 PLANE WAVES 

The equation of motion for a linear anisotropic elastic solid in 
the absence of body forces is given by 


C u - p u = 0 
i j k 1 k , i j i ^ 1 1 


Consider a displacement vector in the following form 


u(x, t ) 


a f 


n.x 'I 

^ J 


( 2 . 21 ) 


( 2 . 22 ) 


in index notation u (x,t) 

i 

where 




n. X 


V 


] 


a = amplitude of displacement vector 
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V = phase velocity 


Substituting {Z2) into (22) yields, 


C n n 
i Jki j 




- P a^. 


f (5) = 0 


t1 

Since function f is arbitrary, f (5) 0, therefore, 







J 


a 

k 


0 


{2.23 ) 


(2.24) 


where F = C n n known as Green-Christof f el Constants. Equation 

ik ijkl j 1 

(2.24) gives a system of three homogeneous linear algebraic equations 
for the unknowns, the scalar v and the vector a . The matrix 
from (2.24) is 


pv^S. 


- r 


i k 


^ a, ^ = <! 0 ^ 
^ J ^ J 


(2.25) 


i.e.. The three eigen values of the phase velocity v and three 
associated eigen vectors for the polarizations can be obtained from 
equation (2.25). Since I is symmetric and real, equation (2.25) 
must have three real eigen values and three linearly independent eigen 
vectors. Hence pv~ satisfies the characteristic equation 


det 


r 


pV-S 


i k 



(2.26) 


The three eigen values of (2.26) gives the possible- wave speeds in 
a given direction of propagation n, and the corresponding eigen vectors 
represent the direction of particle movement. 


2.3.1 Isotropic Solids 

For the case of isotropic solid, the elastic stiffness tensor is 
given by (2.5). So 


r = ( X + p )n n + p n n 5 

1 k is ^ ^ ^ 


(2.27 ) 
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1 , and so 


Since n is a unit vector n n = 

i i 


r-(X+^)nn+u 6 (2.28) 

1 k 1 H i k 

Substituting the value of r in [2,Z6) and simplifying the 

1 k 

characteristic equation reduces to 


- [ X + 2 pv- - ^ ^ = 0 


(2.29) 


Thus, the two possible wave speeds in an isotropic case are given by 


1 



and 


V 

2 



C 

t 


(2.30) 


These wave speeds are independent of n, so the two waves propagate 

with the same speed in all directions. It can be easily show^n that the 

particle motion associated with v and v are parallel and perpendicular 

1 2 

to the direction of wave propagation. 


2.3.2 Transversely Isotropic Solid 

If an orthotropic solid exhibits symmetry with respect to arbitrary 

rotations about one of the axes, then it is called "transversely 

isotropic" solid. The linear stress-strain law for a transversely 

isotropic solid may be given in various forms. Without loss of 

generality, by assuming the symmetric axis of the material as coinciding 

with the X - axis, the elastic stiffness tensor may be represented as 
1 

146J 


23 



i i k 1 


S 1 + . 5 ..+ 6 -S. 5 ^ I 

Ijrni tlikim ikiml Iklmlii i1»1kmi 


^ f' ^ r ^ 

I |5- 5 , 5 s., + + 5. ,5 5. I 

tj\^ilkljm ilmljk jlkliHl jlmlikJ 


+ ^ 6 6 6 6 

it i 1 k 1 m 1 


(2.31) 


The stress - strain relation takes of the form 


~ ^€:..5.. + 2;,j G. + cxf5.,5.,e, , + 5, ,5 .e, 5 .1 

11 kkji til lililkk klmlkmii) 


" h.'l fs. ,5_e, . + 5., 5, ,e. .1+05. ,5 . ,5, .5 f 

tj\^ilklkj jlklkij iljlktml 


km 


(2.32) 


Where X, , a and 0 are the generalized lame’s constants. The 

components of [ r ] can be calculated from (2. 31). The results 

} i 

” u 

are 


= ( X + 2 a + 0 + 4 ;j^ - 2 n^ ) n^ + ^n|+n|) 

^23 = h, n 2 + ( X + 2 ) ng + n^ 

^33 " ^ ^ + 2 ) n2 

X +a + ) n^ n^ 

r, 3 =( \ +a + ) n^ n^ 

^ ^ I ^ 3 (2.33) 

The determinant matrix takes the following formi 
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b n- + b ( l-n~ } -V" 
? 1 B 1 

b n n 
.3 1 7 

b n n 
Pi 1 3 


b n n 

3 1 ? 

b n^+b n^+ b 
B 1 1 ? 4 3 

( b -* b ) n n 
1 4 ? 3 


b n n 
? 1 3 

( b - b )n n 
1 4 3 3 

b n^+b n^+b n^-v^ 
K 1 47 1 


= 0 (2.34) 


where , 


X + 2 /j 

b = L 

1 n 


b = 

7 


X+2a+p+4/Jj - 2fj, 


X + a + /J, 


b = 

■T 


b = 
4 


b = 
s 


(2.35) 


Equation (36) can be factorized resulting in the two equations, 


b n- + b ( 1 - n; ) - V- = 0 

51 4 1 


and 


( b - V- ) + ( b - V-" ) 

B 5 


( 2 . 36a ) 


( b. - 


)n- + (bj 


b )(l-n- 

B 1 


I” ( b,- b ) ( b^ - ) - b2 ] n2 


1 - n- ) = 0 
1 


(2.36b) 


One possible wave speed in the direction n is given by, 


(n) = f n2 + b^ ( 1 - np j2 


(2.37 ) 


The other two wave speeds are the positive roots of equation 
(2.36b). Once the velocities are known, the direction of particle motion 
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associated with each wave can be obtained from (2.30). It can be 
showed that the particle motion associated with the speed given by 
equation (2.37) is always transverse, where as those associated with 
other two speeds are neither longitudinal nor transverse. A detailed 
discussion of the properties of the waves for general direction of 
propagation was given by Hearmon L 4 7 J 


For a typical graphite/epoxy composite containing 60% graphite 
fibers by volume, the Engineering constants are 


= 160.73 GPa 

11 

C =13.92 

22 

GPa 

c = 

IZ 

6.44 GPa 

= 6.92 GPa 

23 

C =C =7. 

55 6 5 

07 GPa 

c = 

44 

3.5 GPa 


In the units of mm/ps the wave velocities for propagation in the 
xz -plane are calculated to be 


V =10 

. 1 

rH 

CM 

II 

> 

II 

> 

11 


Zd 31 

V =2. 

97 

V = 1.49 

31 

33 

A polar 

diagram of 

the wave speeds as functions of propagation 


direction is sketched in fig 2.2. 



Fig. 2.2 Wave Speeds as Functions of Angle of Propagation in 
Graphite - epoxy Composite. 


2.4 Propagation in a half-space 

\vhen an elastic v/ave encounters a boundary between two media, 
energy is reflected and transmitted (refracted) from and across the 
boundary. If the boundary is a free surface, pure reflections will 
occur. What distinguishes elastic wave-boundary interaction from other 
related fields is the phenomenon of mode-conversion. This describes the 
behaviour by which an incident wave, either compression or shear, is 
converted into both types on reflection. 


Rayleigh surface Wave 

In an unbounded elastic media only two types of waves can be 
propagated. However, when there is a boundary, as in a half space 
problem, a third type of wave exists w^hose effects are confined closely 
to the surface. These are called Rayleigh surface waves and their 
effects decrease rapidly w^ith depth and their velocity of propagation 
is smaller than that body waves. 

By investigating the propagation in a half space the frequency 
equation for surface wave may be obtained as, 

s! r (£_)6_ 8( — )"+ (24 = 0 (2.38) 

2 C _ C ^ 

C ^ 2 2 ^ 

? 


Where c is the surface wave velocity and k- 

Out of the three roots of the equation, it has been shown that for 
all real media (0 <p< 0.5) there is only one real root. The resulting 
surface wave propagating with the velocity established from the above 
equation is the Rayleigh wave. An approximate expression has been 
developed for Rayleigh wave velocity, which is slightly less than shear 
wave velocity, given by, 


( 0.87 + 1.12 V ) 
( 1 + p ) 


(2.39) 
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The particle motion can be shown to be elliptical in nature and 
retrograde with respect to the direction of propagation (i,e, it is 
counter-clockwise for a wave travelling to the right). The amplitude 
decreases exponentially away from the surface. Rayleigh waves are 
essentially tw'o dimensional. Hence, energy associated with these waves 
will not disperse as rapidly as that of three dimensional waves of 
dilatation and rotation. 

2.5 Waves in Layered Media 

In many situations waves originate and propagate in a media having 
a layered structure, where interface between dissimilar materials exist. 
An important example of such a case is composite laminates. The 
impedance mismatch (the product of density and wave velocity) between 
the adjacent layers is responsible for the reflection, refraction and 
transmission of the waves. Interface waves may also generated. It has 
been shown that the existence of such a wave required that the 
shear-wave velocities of the two media has to be nearly the same. These 
interface waves are usually referred to Stonley w'aves in the literature. 

Through composite materials are inherently inhomogeneous and 
anisotropic, when the scattering by' defects which are large compared to 
the fiber dimension is considered, it may be assumed to be homogeneously 
orthotropic. In this case, the wavelen'gths involved are much larger than 
the fiber dimension. 

2.6 SCATTERING OF ULTRASONIC WAVES 

In this section the formal aspects- of wave scattering by a single 
flaw are discussed. These aspects are independent of the specification of 
the inhomogeneity/flaw. 

The scattering problem refers to a solution of the field equations 
in the regions interior or exterior to the volume V of the scatterer. 
The scatterer is excited by a source exterior to V. The field produced 
by the source in the absence of the scatterer is called the incident 
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lieid The dillereiice between the total field outside V in the 

presence ol the scatterer and is called the scattered field 
- All three fields, and separately satisfy the field 

equations in their respective regions of validity. 

Scattering may be viewed as a two step process. According to 
Huygens principle, the incident field excites sources on the surface S 
of the scatterer, which in turn radiate to the exterior and the interior 
producing the scattered field and the interior field respectively. Since 
the sources of the scattered field are confined to S, physical 
considerations require that u' must be finite and decay in a particular 
manner at distances far from S. This requirement is called the 
Sommerfeld radiation condition. A solution of the field equations, that 
simultaneously satisfies the boundary conditions on S as well as the 
radiation conditions at infinity, formally defines the scattered field. 
If the known sources are inside V or distributed on S, then the 
corresponding problem is called the radiation problem. 

The basic scattering problem is shown in fig 2.3. In a material 

having elastic stiffness constants C and density p, there exists a 

i j k 1 

flaw defining a region of space R and having a surface S. The material 

} j 

inside the flaw has stiffness constants C and density p w^hich, in 

i i k 1 

general, differ from those of the host material. For many cases of 

$ ) 

interest, C = p = 0 (void). In ultrasonic flaw detection 

i ikl 
^ ^ - 

experiments power generated by a piezo-electric transducer through a 
longitudinal/shear wave is directed towards the flaw. The presence of 
the flaw, however, prevents the incident power from propagating 
unhindered, and flaw selectively scatters part of the incident power 
into various directions. The transfer of energy takes place at least in 
both shear and compressional modes. In addition it was found that part 
of the incident wave energy striking a stress-free cavity is converted 
into a surface or boundary wave which " creeps " around the cavity 
radiating energy in all directions (like a light-house beam). The 
measurement of the scattered power into a given direction, relative to 
the incident power is of practical interest. Such a measurement may 
serve to differentiate between flaw geometries, and hence, would be an 
useful ultrasonic NDE measurement. 
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The physical quantity that measures the scattering effectiveness is 
the " cross section " of the flaw. This is related to, but not equal to, 
the effective geometrical cross-sectional area of the flaw. 


An intensity vector 1 associated with a stress field a and a 

i 11 

t " “ 

velocity field u. may be defined as, 


I = a u (2.40) 

t ill 


The term intensity means the rate at which transfer of energy per unit 
area normal to the direction of propagation, i . e. , 


I = I n (2.41) 

i i 

Where, n. is the unit normal to the direction of wave propagation. 
For a plane wave incidenting along the z-axis, the incident intensity is 
calculated by considering the rate per unit area in which energy is 
carried across an infinite plane perpendicular to the z-axis. Thus, 

a^u^ (2.42) 


The power flux or rate at which energy is carried across a surface A 
with unit normal n is given by, 

i 

. P = j' n. dA 

=rnaudA (2.43) 

.1 i i i i 


For a given angular frequency, corresponding to a period T, the total 
cross-section P(o)) is the ratio of the average power flux scattered into 
all directions to the average intensity of the incident fields. 


P(w) 


< P- > 

< I > 


(2.44) 
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Where the angular brackets denote the time averaging, 


1 ’’’ 

< f(t) > = -i- r f(t) dt 

'^0 (2.45) 

For a given frequency, a measure of the fraction of incident power 
scattered into a particular direction is knowm as " differential cross 
section " and denoted by ■ where, dQ is the differential element 

of solid angle. 


The average power flux, relative to the average total incident 
intensity directed across an element of area dA normal to n is given 

i 

by, 


dPIo)) 


< n a u- > dA 
- - j - 


< l” > 


At a distance r from the origin, dA = r“do. Thus, 


dP(w) 

dO 


.2 • s . 

< r n <7 u > 
i i j J 


< 1 ° > 


(2.46) 


(2.47) 


In more precise terms, 


JT-., V 1 • < r n a. . u > 

dP((o) _ lim > 1 J J (2.48) 

dQ r-» on ^ jO ^ 

The limiting process is necessary since differential cross section 
measurements are made at distances very far from the scatterer (the 

flaw). It has been found that as r-^ oc. , a® . and u® are of the order of 

i J J 

1 . Therefore, as r cc , equation (2.48) represents that portion of 

r 

the scattered power in do which is independent of the distance from the 
flaw . 
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Gubernatis et ai . [48 J have showed that the general expression lor 
differential scattering cross section for a plane wave incidenting on an 
arbitrary shaped flaw in an isotropic medium as 


dp{(x)) a (X+2^) Bn 

do « (X+2,j) |a .(2+ (b j.^ 

where a and p are the square of the 
shear wave respectively, |A j and |B. j 
ja^l and j | are the amplitudes of 
waves . 


wave numbers of longitudinal and 
are the scarred wave amplitudes, 
incident longitudinal and shear 
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CHAPTER III 


FORMULATION AND COMPUTER IMPLEMENTATION 


This chapter deals with the finite element formulation of the 
governing wave equation which is a hyperbolic partial differential 
equation. Spatial and temporal discretization problems, accuracy and 
stability requirements were discussed. Transducer effects were taken 
into account to some extent. Methodology of differential cross-section 
measurement is given. A brief account of computer implementation of the 
model is discussed. 


3.1 FINITE ELEMENT FORMULATION 

3.1.1 FORMULATION 


The governing hyperbolic partial differential equation for a 
linear homogeneous and anisotropic elastic solid neglecting body forces 
is, 


C 


u 

i jkl k 


1 j 


. “ P 


a'u. 


at' 


(3.1) 


The total potential of a linear elastic body using D* Alembert 
principle is given by 




1 

' a"u. 

u 

i 5 J 

'J ’.r 

a . dV - 
1 J 

u t dS + 
i i 

C* *■ 

u p ■ dV + 

at" J 


1 

3t 


(3.2) 


where the first term is the strain energy, second term is the work 

done by the surface traction t and the third term is the work done by 

i 

the inertia forces and the last term is energy dissipation due to the 
internal friction with g being coefficient of damping. 


The finite element solution involves the discretization of the 
domain into a number of elements, approximation of the field values 
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interior to the elements in terms of its nodal values through the shape 
functions ol the chosen element and the determination of the nodal 
values through the minimization of the above energy functional. Invoking 
the first variation of n, i.e, 611 = 0, the following equation is 
obtained 





5u a dV - 

i > j i j 

5u t dS + 

6u pu dV + 

i i 

1 i 


. g 

V 


5u . ^u , dV = 0 


(3.3) 


For the field variable u in eqn.(3.3) it requires continous 

I 

element. Approximating the field values u with in an element in terms of 

I 

the nodal values (U) as 




NU( t) 




( np ) 


where 




is the set of shape functions. 


The strain-displacement now, is represented as 


(3.4) 


(x,y, 2 ,t ) J. = 


1=1^ B^-^^(x,y,z)j|u{t) 


( np ) 


where 


[B^] 


is the strain-displacement matrix. 


(3.5) 


The constitutive equation is 


^a^®^(x,y,z,t) 1. = r D^®^l^e^®Mx,y,z,t) j. 

L ; L n J 


(3.6) 


where is the material property matrix and depends on the problem 

type whether it is plane stress or plane strain. 


Combining equations (3.3) - (3.6) and several algebraic 

manipulations gives the equation of motion in the matrix form as 
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(3.7) 


r -] r ^ 

' V / ' I I 


M 1 -i ^ M I c I ^ }. + I K 1 .{ u 


I at^ j 


I, at j 


I J 


r ^ 

< F I' 
I J 


where , 


r M 1 
I « 




I ® r r 

M = I p N 


-,T 


1- J *—1.-1 


Jv 


r 1 

N dV 


U -J t- J 


ri vri®r ri^rn 

1 C I = I 5 [ N I I N ] dV 

-LJ Jy 


C 

I i 

U J 


I K I 


vf 1® I’ r I’^r ir i 

= ; I K I = I I B I I D I j B I dV 


L J L. -i L -r 


<! F }. = ) ^ F ;. = I N .{ t J. dS 

L J “ L J Js L J L J 


(3.8) 


The matrices [Mj, [C], [Kj and {F} are well known as mass matrix, 
damping matrix, stiffness matrix and forcing vector, respectively. 
Mass, damping and stiffness matrices are sparse, symmetric and 
positive-definite. Since the matrices are very large and sparse, so 
they were condensed efficiently by skyline storage algorithm. 
SYni-th et al [31] have examined the accuracy of finite element 
solutions and found a ringing phenomenon which dies away only 
gradually. These oscilltions can be suppressed efficiently in the 
presence of damping matrix. But it is a tedious job to find out ^ 
distribution. It was observed that in wave propagation problems the 
frequency of ringing is related to the size of the elements used in 
the FE mesh, and found to be predominant at frequencies of about4-6 
nodes per wavelength. But perfectly elastic formulation is 
incorporated in this analysis by taking 10 nodes per wavelength. 
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3.1.2 Spatial Discretization 


For wave propagation problems low-order elements were recommended 
this reason linear quadrilateral isoparametric elements 
have been chosen for spatial discretisation, a typical element is 
sketched in fig (3.1;. 


V 



(c) 

Fig. 3.1 (a) Natural Coordinate and Rectangular Parent 
Element. 

(b) Isoparametric Element (Mapped on Cartesian 
Coordinate; Quadrilateral Elemet.) 

(c) Shape Function. 
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3.1.3 TIME DISCRETIZATION 


The system of linear, second order, and ordinary differential 
equations with constant coefficients (eqn. 3.7) has to be integrated 
in the time interval {0,T) by discretising the time variable as t = n 

n 

/it, 0:^ n £ N where At = T/N. 

NEWMARKS METHOD 

It was reported that Newmark’s algorithm is more accurate than 
cental difference scheme and therefore it is adequate to understand the 
physical phenomenon of wave propagation and to study arrival times and 
amplitude effects [35]. The Newmark’s constant-average acceleration 
scheme makes the following assumptions, 


t+At 


au 

at 


^au 


+ f (1-5) 

V. 


at- 


+ 5 




At 


at’ 


(3.9) 


t+Aty _ t-y 


■au 




At + I (0.5-a) 


at“ 


+ a 


t+At^2 
o 


at 




At 


where a and 5 are parameters that can be determined to obtain 
integration accuracy and stability. Newmark [49] originally proposed an 
unconditionally stable scheme known as constant average acceleration 
method, in which a = 0.25 and 5 = 0.5 


In addition to the above 
displacements , velocities and 
at time t+At are considered: 

M , K 

at“ 


difference equations, for solution of 
accelerations the equilibrium equations 

= (3.10) 


Using equations (3.9) in (3.10) we obtain the following equation 


for 
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t.Aty = _|tf 




(3.11) 

It is seen that this is an implicit scheme and requires the 
factorization of the stiffness matrix. 


3.1.4 ACCURACY AND STABILITY CONSIDERATIONS 

Numerical modelling of elastic wave equation introduces several 
errors due to time and space discretization. These errors can be 
minimized by the careful selection of spatial grid size and time step. 
At least eight nodes per shortest wave length are chosen to minimise the 
error with the spatial resolution. In order to represent all the 
frequencies to pass through the mesh uniform low order elements are 
chosen . 

The time step size should be such that the wave travels less than 
the spatial grid size in one time step. In other words, the energy 
should not propagate right across one mesh volume in one step. 

The stability conditions are related to the wave velocities by, 


At = 


(3.12) 


where h is the minimum spatial resolution and c is the maximum 

m 

propagation velocity, i.e, that of dilation wave velocity for isotropic 
solids and maximum phase velocity of quasi-compression waves for 
anisotropic solids. Another criteria suggested are 127] . 


At ^ 


/ ^ A. ^2 ^ 0:5 

( c + c ) 

1 ?. 


for isptropic solids, 


(3.13a) 


where c and c are the compression and shear wave velocities. 
1 ? 


At £ — =T 77 : for anisotropic solids. 

1.732 c 


3 . 13b) 


39 



whei’e is the v'elocity of quasi-compression waves. 


3.2 BOUNDARIES OF THE DOMAIN 


The physical region modelled is very small compared with regions 
considered in any experiment. It was observed that in numerical 
modelling of wave propagation problems the presence of artificial 
boundaries introduces spurious reflections which contaminate the 
solution. The problem in this present work has been overcome by 
constructing a domain of sufficient size by properly locating the flaw 
so that the required solution is obtained before the reflections from 
the external boundaries arrive. 

In the literature nonreflecting, absorbing and transmitting 
boundaries were proposedL 25 , 50 , 51 J , so that very little reflections 
takes place. These techniques involve the superposition of solutions 
with different boundary conditions and thus much more costly than a 
single solution with respect to computational time and cost. 

3.3 MODELLING OF INPUT PULSE 

A main step towards the successful modelling of elastic fields 
associated with NDE is the generation of realistic pulsed ultrasonic 
waves through the use of piezo-electric transducers. Neglecting the 
internal transducer mechanisms, the necessary transient excitation 
pulse can be specified by either force or displacement variations 
across the aperture of the transducer. Of these two possible choices, 
force type boundary conditions are chosen in this present work. These 
force excitations will have an amplitude and polarizations imposed on 
them to represent the transducer pulse. 

\ 

The forcing vector chosen is specified as shown in Fig (3.2) and 
mathematically it can be represented as, 5 
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F(x,y,t) = p a 6(x) w(y) f(t-0{y)) 


(3.14) 


where 

p is the polarisation associated with F 
a is the amplitude of forcing vector 

6(x) ensures that F is applied on the surface of the 
solid only. 

w(y) takes into account the possible non-uniform 
distribution of the force across the aperture. In other words, it acts 

as a weighting function a uniform distribution w(y)=l 

for — ^ < y < L- is presently assumed. 


f(t-0(y)) stimulates a realistic transducer signal that is to be 
determined experimentally by measuring the received transducer voltage 
signal in a pulse-echo experiment. The signal is subsequently digitized 
to be used as an excitation. A good analytical approximation is given 


by[ 52 J , 


r 

I 

f(t) = \ 

I 

I 


j" 1 - cos(a)t)"| cos(3£jt} 

V j 

0.0 


t i 1 /jS 

t > 1 ;.jS 


(3.15) 


where, g) = 2jr (transducer frequency). This approximation models the 
transducer with a harmonic time dependency modified by a cosine-shaped 
weight. The possible phase delay 0(y) for beam steering purposes is set 
to zero. 


The boundaries away from the transducer aperture are set to 
traction-free conditions i.e., <7,.n. = 0, 

I 1 I 

The initial conditions for starting the iterations used are, 


U(t =0) = 0 

U(t=-At) = 0 


au 

at 


(t=0) 


0 


(3.16) 
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W(y) 


Y 



Fig. 3.2 Specification of force vector on boundary 
to simulate transducer pulse 

3.4 DETERMINATION OF SCATTERING PARAMETERS 

Estimates of differential scattering cross sections can be easily 
obtained with this model. To obtain these cross sections, the 
calculations are first performed with the flaw present for the given 
source position and its polarisations. Reference calculations are then 
, 'phese are the same as the original calculations in terms of 
mesh size used in the domain, time step, source location and duration of 
pulse application, but differ only in the absence of the defect. The 
vector displacements from the appropriate reference calculation are then 
subtracted from each original calculation at every mesh point and at all 
time steps. The resulting displacements are entirely due to the 
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scattering effects of the defect. 


Mathematically , 

i 9 i _i • 1 

u " ■ „ displacement field obtained with the defect 

' » incident field (obtained when there is no 

defect ) 

u”"- = scattered field 
= . u'*’ 

The energy in the various pulses is proportional to the square 
of the displacements and so, it is a simplified measure of the energy 
scattered. So, the expression for the energy scattered is given by 
L28j , 


t 

p'2 

® = (■ + vO dt 

■“ t '' ^ 

1 

where u , v are the displacements at a detector. The integral is 
performed over the whole time trace i.e., from 0 to the A-scan limit. 
The integral here is performed using Simpson’s rule. The value of 
displacements u, v at any point in an element are interpolated from tlie 
nodal values. 


3 . 5 COMPUTER IMPLEMENTATION 


The two-dimensional finite element code has been implemented on 
CONVEX C220 to model pulsed ultrasonic wave propagation in linear 
elastic solids with a scale of inhomogeneity. The block structure of the 
program consists of three parts : a preprocessor, a processor and a 
postprocessoi' . 


The preprocessor generates the required data for the 
discretization of the domain. It consists of co-ordinates of nodes and 
the nodal connectivity of quadrilateral elements used for the 
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discretization* The material properties, problem type, transducer 
characteristics are also included in it, A graphics program was 
developed for drawing the finite element mesh and it also checks the 
connectivity of elements. 


The processor is the main body of the code and models the wave 
propagation. It will predict the vector displacements at all nodes in 
the domain at all time steps using Newmark*s algorithm* Separate 
routines were incorporated to calculate the strains and stresses in 
the elements. These can be used to find out the strain in the elements 
at different time instants. Similarly kinetic energy can also be 
calculated . 

The postprocessor is a starbase graphics program. It visualizes the 
complete scattered wave field at any desired instant. A scan data can 
also be plotted. Differential scattered cross-sections can be drawm with 
graphics code. 
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CHAPTER IV 


RESULTS AND DISCUSSION 


This chapter presents the details of various simulations 
performed on homogeneous and composite materials with the Finite 
Element code being .described in the earlier chapters. These results 
are sho-wn in the form of A-scan plots, differential scattering 
cross-sections and complete wave field visualizations at different 
instants {snapshots) and thus giving ultrasonic testing. 

4 . 1 RESULTS FOR ISOTROPIC SOLIDS 

Three cases have been studied for the homogeneous and isotropic 
medium. For each case both the longitudinal pulse and shear pulse 
interactions with the flaw are studied. Aluminium is chosen as the 
material. In all the finite element calculations only one half of the 
domain has been discretized taking the advantage of symmetry. 

Four nodded quadrilateral elements for spatial discretization and 
Newmark’s algorithm for time marching were chosen in the present 
analysis [38]. Numerical experiments have been conducted to select the 
optimum spatial discretization and time step values. For selecting the 
maximum spatial discretization value the pulse arrival time at a 
particular location is compared with the theoretically computed value. 
The time step value has been obtained from the equation 3.13. 

The Finite Element main processor output includes the x and y 
displacements at all nodes for each time step. These results are 
presented as "snapshots". A snapshot represents the displacement field 
( either x or y ) at a particular time step, which is plotted normal 
to the plane of the two - dimensional domain. Incidently this type of 
presentation shows the mode-conversion, diffraction and interface 
effects which are also observed in photoelastic experiments. 

The scattered amplitudes of the displacements are obtained by 
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subtracting the incident field, which were obtained from a defect 
free plate, from the total field. In order to minimise the 
computational time, uniform mesh has been chosen for the defect free 
plate. The displacements at the corresponding nodes of the defective 
plate are calculated using interpolations. But some spurious 
displacement field is noticed in the snapshots which were taken before 
the start of scattering time. This important observation may be 
attributed to the dispersion of wave when it propagates through the 
finite element mesh. This phenomenon can be explained in the following 
physical terms: 

Any transducer pulse signal contains all frequencies whereas the 
finite element mesh has a limit onf requencies it can describe 1.53J. 
Thus a pulse tends to broaden as it propagates, if there is an 
artificial boundary, caused by the expanding cell, there is reflection 
of the dispersing wavefront, which leads to noticeable effects. (An 
element which is incapable of describing a particular frequency may 
also act as an artificial boundary). In the absence of these 
reflections, dispersion leads to a stable waveform which is formed 
from the given input wave profile after a small number of time steps. 
This new waveform propagates without further detectable change in 
shape. In order to minimise this peculiar effect of wave propagation, 
the domain discretization of defect free plate is closely maintained 
with that of the plate with a defect. 

In the present analysis the differential scattering 
cross-sections are plotted as the angular distribution of scattered 
component of wave amplitude. For this reason these plots can be better 
termed as a polar plot of scattered amplitudes. The angular 
distribution of scattered component of wave amplitudes are determined 
by observing the displacements at three to four mesh spacesaway from 
the maximum size of the flaw. Final traces are obtained by integrating 
the u and v components of scattered displacements over the entire 
A-scan period, so as to cover the first few scattered signals. 
Obviously, waves will propagate back and forth along the flaw, so that 
some scattered energy will arrive later; it is in this respect that 
the present analysis differs from the analytical calculations. 
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The elastic constants of Aluminium are taken as : 

E = 70 GPa u = 0.33 p = 2275 Kg/m' 

The calculated compressional , shear and Rayleigh wave velocities 
in A1 are : 

c = 6750 m/sec c = 3375 m/sec c = 3150 m/sec 

1 ? R 

TABLE 4.1 


FE simulation details for various defects in isotropic solids 


Problem 

No. of 

No. of 

At 

F 

Flaw Size 


elements 

nodes 

/jsec 

M Hz 

mm 

1 . Defect 






free 

7500 

7676 

0.04 

0.75 

— 

2 . Surface 






Crack 

7500 

7691 

0.04 

0.75 

15 

3 WEDGE CRACK 

- 





a. With wedge 

6656 

6839 

0.04 

1.0 

15 

b. Defect 






free 

6656 

6825 

0.04 

1 

1.0 

1 

1 - - - ... 

— 

4 . Drilled 






Hole 

1 

i 





a. With Hole 

5722 

5887 

O 

• 

o 

1.0 

2.4 dia 

b . Defect 

5632 

1 

I 

5785 

0.04 

1.0 


free 






Care has been taken to use about twelve nodes per longitudinal 
wavelength. The domain modelled is of 30 mm x 80 mm for the defect 
free and surface crack modelling, 31.2 mm x 44.4 mm in the case of 
awedge crack and 35 mm x 50 mm in the case of a hole. In each case of 
the defect, the pulse is polarised along x-axis and then along y-axis 
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to resemble the longitudinal and shear pulses. 


4.1.1 Propagation in defect-free isotropic solid 

The geometry modelled is as shown in fig. 4.1 along with a 
typical finite element mesh. The transducer pulse width is assumed to 
be 4.0 mm and is uniformly distributed over the surface nodes situated 
above and below the symmetric portion. 

Figs 4.2 and 4.3 shows the snapshots of X and Y displacement 
fields at different instants of time after the application of pulse 
(Ifjsec). At 1.8 fjs , which is just 0.8 ps after the application of 
pulse, and the details were not clear. The next three plots (3.6, 5.0 
&. 6 . Opsec ) clearly show the fast and radially propagating compression 
wave front (Pressure wave, P) and the compressed or expanded parts are 
also seen. Outside the source no deformation exists. A very strong 
shearing deformation is observed {rotational wave, S) between the 
deformed area of the longitudinal waves and the undeformed area. These 
transverse waves propagate obliquely. Also a Rayleigh (Rl) wave 
propagating along the surface with a velocity slightly less than that 
of the shear wave is seen. At the instant of 6.0 ^sec the compression 
pulse has reached the back wall and is reflected with a negative 
reflection coefficient and is travelling towards the front side of the 
block . 


In fact, real transducers produce more than one type of waves 
that they are designed to give, e.g. compression wave transducers also 
produce shear waves. So, the use of the waves generated by all the 
components of the incident wave field provides more information for 
defect characterization [26] . 

4.1.2 Rayleigh wave interactions with a surface crack 

Surface breaking cracks (such as fatigue crack, debris )are the 
most commonly encountered flaws on the structural components. In the 
present analysis this crack is modelled as a narrow cut of essentially 
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zero width in the mesh. A schematic representation of the domain 
considered and the typical finite element mesh are as shown in the 
fig. 4.4. Two surface cracks assumed to be present in the model to 
make the domain symmetric. 

It is observed from the previous result that a compression wave 
produces a strong surf ace/Kayleigh wave. Figs 4.5 and 4.6 show the X 
and \ displacement fields, when a wave polarised in X— direction 
interacts with the surface crack. Strong scattering of Rayleigh wave 
by the surface crack can be seen from the shots of 5.0 /Jsec and 6.0 
psec. Only a small part of the incident Rayleigh wave is transmitted 
beyond the surface crack. 

Figs 4.7 and 4.8 show the waves generated entirely due to the 
scattering of Rayleigh wave by the surface crack. (These are obtained 
by subtracting the incident displacement field from the total field). 
As the incident Rayleigh wave reaches the crack corner, a part of its 
energy is reflected back to the transceiver (in pulse-echo mode) and 
the rest travels down the crack face as a Rayleigh wave. In addition 
there is also mode-conversion to compression and shear energy from the 
crack corner. The Rayleigh wave travelling down the crack face 
mode-converts itself to a shear wave at the crack tip which propagates 
into the material. The shot at 6.0 p sec clearly visualises this and 
the radial distribution of mode-converted compression from the surface 
crack corner. This mode-conversion is in agreement with the 
experimental observation of Bond [26j. A schematic diagram of this 
Ultrasonic NDT is sketched in fig. 4.9. 

A polar plot of scattered amplitudes for this case is shown in 
fig 4.10. The amplitude is seen to be maximum in the forward 
scattering direction rather than in the backward direction (180"). The 
lobe structure indicates that the apparent area, blocking the incident 
wave, in forward-scatter is much more than the back-scatter. 

4.1.3 Compression wave interaction with a surface crack 

Figs 4.11 and 4.12 show the X and Y displacements when a wave 
polarised in a direction perpendicular to the propagation direction 
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{ SV wave ) interacts with the surface crack. The snapshot at 3.6 p sec 
clearly shows the compression wave interaction with the surface crack. 
The strongly reflected Rayleigh wave and shear waves have in turn 
interacted with the incoming Rayleigh and shear waves respectively. 
The final shot shows the reflection of the pressure wave from the back 
wall. In the Y- displacement, the 6.0 fj sec shot shows the reflected 
compressive wave and the generated Rayleigh wave on the back wall. 

The scattered X and Y displacements for this case are as shown 
in figs 4.13 and 4.14. Because of the interaction of compression wave 
with the surface crack, the strongly reflected compression, shear 
waves which are propagating into the material and the reflected 
Rayleigh waves which are propagating along the surface are clearly 
revealed in these snapshots. From the final snapshot, one can clearly 
see, that the angular distribution of the scattered shear wave, which 
originates at the crack tip, and also the curvature of the pressure 
wave front. This observation is in agreement with the photo elastic 
experiments [53J and a line diagram of the scattered wave fronts is 
shown in fig 4.15. 

The reflected Rayleigh wave seems to be travelling outwards from 
the sides of the surface crack. In a practical situation, if two 
receiving Rayleigh wave transducers are placed on either side of the 
surface defect, the relative amplitude and the shape of the received 
Rayleigh pulses at the two receivers contain information about the 
orientation and size of the surface crack 126]. A schematic 
representation of this is as shown in fig. 4.16. 

Fig 4.17 shows a polar plot of scattered amplitudes when a 
compression wave interacts with the surface crack. The curve seems to 
be similar in ' the forward and bcak-scatter except for their 
amplitudes. There is an enormous decrease in the scattered energy in 
the 180" direction compared to the O" direction. As expected a 
relatively low amount of energy is scattered in 90'’ direction. The 
energy is likely to be maximum in the 65* direction in both forward 
and back-scatter. 

4.1.4 Compression wave interaction with a hole 
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Cylindrical defects like side drilled holes, voids are another 
class of structural encountered flaws. Fig 4.18 shows a typical 
finite element mesh generated to incorporate this type of defects. The 
hole is ot 2.4 mm in diameter. The transducer pulse is of 3.2 mm in 
width . 


Figs 4.19 and 4.20 show the X and Y displacement fields when a 
pulse, which is polarised along X-axis is applied, interacts with the 
hole. The differentiation of pressure and shear pulses is very clear. 
The complex interaction of Compression pulse with the hole is observed 
in the snapshot at 3.8 psec . The slightly reflected compression pulse 
is also seen here. The snapshot taken at 5.0 /jsec clearly visualises 
the incident wave and scattered wave fields. 


Figs 4.21 and 4.22 shows the scattered X- and Y- displacement 
fields. There is no scattered field till 2.8 ^sec shot as there is no 
interaction of the wave with the hole. A high and complex scatter 
field is observed surrounding the hole at 3.8 psec shot. The 
scattered pressure pulse, which is propagating at a higher speed and 
radially distributed from the center of the hole, is followed by the 
mode converted shear pulse. The distinction between scattered 
compression and shear wave fronts is predominant in scattered Y- 
displacement field. The mode converted Rayleigh wave is not clearly 
seen in these plots. 


The angular distribution of scattered component of w^ave amplitude 
is as plotted in fig 4.23. There is a sharp lobe structure in both the 
forward-scatter and back-scatter directions. The main contribution to 
this plot seems to be due to the scattered compression wave as it is 
also radially distributed. 


Fig 4.24 shows the A-scan display at a point before the hole and 
45- with the line of propagation. The compression pulse expected 
arrival time is 2.5 /jsec. The reflected pressure pulse, the mode 
converted shear pulse and the direct shear pulses can also be seen in 
the figure. 
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CENTRAL L??RARV 
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A. 




When the wavelength of the elastic wave is approximately equal to 
the diameter of the hole, the incident pulse would be expected to 
excite normal modes of vibration in the target. In order to study this 
a hole of diameter equal to 6.4 mm in a domain of 45mm x 58 mm is 
modelled. Figs. 4.25 and 4.26 show the X- and Y- displacement fields 
and Figs. 4.27 and 4.28 show the scattered X- and Y'- displacements. 
When compared with the earlier case the scattered displacements are 
showing a wave ( surf ace ) which is creeping around the circumference of 
the hole and radiating energy in all directions. These creeping waves 
at an eigen-frequency match phases during repeated excursions round 
the hole circumference, and by constructive interference build up to 
resonance [26] . The effect of this is felt in the differential 
scattering cross-section shown in fig. 4.28, where two distinct lobes 
are seen. The scattered energy seems to be maximum along 55^ direction 
in forward-scatter. 

4.1.5 Scattering of SV wave with a hole 

To study the shear wave interaction with the hole, the transducer 
input is polarised in Y- direction and the wave is set to propagate 
along X- direction. 

The scattered X- and Y- displacements are shown in figs 4.30 and 
4.31 which are magnified greatly. The scattered field at 3.8 fjsec is 
due to the head wave interaction with the hole. The snapshot taken at 
5.0 /Jsec clearly reveals the radial distribution of mode-converted 
compression pulse and scatted shear pulse. Here also the creeping wave 
distribution is not clearly seen as the surface area available for its 
motion is less than that of the wave length. 

Fig 4.32 shows the angular distribution- of the scattered 
amplitudes .'Here three sharp peaks are observed, two lobes in the 
forward— scatter and one in the back— scatter areas. These peaks may be 
due to the smallsize of the hole when compare with the incident pulse 
wavelength . 
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4.1.6 Compression wave interaction with a wedge crack 

A typical finite element mesh to model a wedge crack is shown 
in the fig. 4.33. The FE mesh is fine near the crack surface and is 
uniform away from the crack tip. The transducer pulse width is equal 
to the wedge length i.e., 15mm. 

Figs 4.34 and 4.35 show the visualisation of X and Y 
displacement fields respectively when the compression pulse interacts 
with the crack surface. At 3.0 psec the compression wave front has 
just reached the crack surface. The snapshots followed by it (i.e at 
3.8 fjsec and 4.8 /jsec) clearly shows the scattered, incident and 
transmitted waves. The edge diffraction of the incident wave can also 
be seen. 

The scattered field is visualised in figs 4.36 and 4.37 from 
the snapshots of 3.0 fJsec and onwards. Since the magnification factor 
used is quite large, so the distinction between reflected pressure 
pulse and mode-converted shear pulse is not very clear in fig 4.30. 
But the strong reflection of compression pulse may be sensed by 
transmitting transducer itself in the pulse-echo mode of testing. The 
distribution of scattered fields is very clear in Y displacement 
field. These scattered waves are sometimes called the edge waves as 
they emanate from the crack tip. The mode-converted Rayleigh wave 
seems to be present surrounding the wedge surface. But the available 
surface for its propagation is negligible , so it may again convert 
itself into the body pulses. 

The angular distribution of scattered amplitudes is sketched in 
fig 4.38. Here three different plots are plotted, each one is computed 
at different radial distance. There is a complete similarity in these 
curve shapes, but the change is only in their amplitudes. As we go 
far away from the crack position, the amplitudes are decreased as is 
generally expected. There is no scattering at all in the 90' 
direction. The amplitudes are maximum at around 50® in both 
back-scatter and forward-scatter directions. 
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Fig. 4.39 shows the difference between the A-scan plots taken on 
the symmetric surface with and without the wedge crack. The direct and 
reflected pulses are noted in the figure. The input pulse amplitude 
is gradually decreasing in the absence of wedge crack at the point 
considered. The effect of input pulse force distribution is clearly 
seen in these plots. 

4.1.6 Scattering of Shear waves by a wedge crack 

In order to study the shear pulse interaction with the wedge 
crack a pulse polarised in Y- direction is applied. The radial 
distribution of mode-converted shear pulse and reflected shear pulse 
can be observed in the scattered displacement fields shown in figs 
4.40 and 4.41. The scattered field is more or less similar to the 
earlier case, except for the difference in the amplitude. 

The angular distribution of scattered amplitudes are plotted as 
a polar plot in fig 4.42. There are two lobes in both forward and 
back-scatter plots. Very less scatter amplitude can be observed in 
both 0*^ and ISO^directions . One can observe that a little amount of 
energy is scattered in the 90^ direction. The area blocked by the 
crack in forward-scatter is much higher than that of back-scatter 
direction . 

4.2 RESULTS FOR COMPOSITE MATERIAL 

Composite materials are anisotropic in nature and in the present 
analysis it is assumed that itis homogeneous and orthotropic, for the 
purpose of studying the wave interaction with flaws which are much 
larger than fiber size. The composite material chosen for the analysis 
is Graphite/Epoxy, the material properties of which are given below. 
Symmetry has been incorporated in all the scattering geometries 
modelled to simplify the finite element discretization. 

The properties of Graphite/Epoxy are, 

E = 181.0 GPa. , E = 10.3 GPa. , G = 7.17 GPa. 

= 0.28 = 0.42 p = 1580.0 kg/m^ 
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The corresponding orthotropic elastic constants for 2-D analysis 


can be calculated to be, 


C =: 181.8 

11 

GPa. , 


c 

12 

= 10 

.34 GPa 


C = 10.34 

zz 

GPa. , 


C 

44 

= 3. 

63 GPa 


C =7.17 

66 

GPa. , 


C 

23 

= 3. 

09 GPa 


C = C = 

16 26 

0.0 GPa 






Due to the 

anisotropy 

of 

the material 

the propagation 

velocities 

for each wave are different 

in different 

directions. The 

analytical 

calculations of 

the phase 

velocities 

f or 

Graphite/Epoxy laminate in 


different directions can be calculated using Green-Christof f el 
equations (Eqn 2.25). Some salient values are tabulated in Table 4.2. 
The calculated values of maximum and minimum phase velocities and the 
corresponding wavelengths of quasi-longitudinal ( QL ) wave are given 
below. The QL wave has a maximum velocity at 90", i.e., along the 
direction of lowest stiffness and minimum at 5^. 


C = 10700 m/sec . , 

m y q 1 


X =10.7 icm. , 

mxn 1 
* - ** *1 ~ 


C = 2468 m/sec 

mnti 1 


X =2.5 mm 

mxn 1 
.... _ 


The corresponding values for quasi-transverse (QT) wave are 
calculated as, 


C =2698 m/sec . , 

mxn t 


X = 2 . 7 mm 

m X n t. 


The maximum is at an angle of 


C = 2130 m/sec 
mnqt 


X =2.13 mm 

m n n 1 


12' and minimum at 0^ and 90'. 
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TABLE 4.2 


Angle : deg 





C :m/sec 

nt 

0.0 

2559 

2130 

50.0 

8337 

2284 

10.0 

1 

2910 

2427 

60.0 • 

9352 

2175 

o 

o 

4221 

2440 

70.0 

10106 

2175 

30.0 

5695 

2400 

80.0 

10570 

2142 

40.0 

1 

7098 

2345 

1 

90.0 

1 1 

10726 

2130 

1 1 


The shortest wavelength was found to be 2.15 mm. for the frequency of 
1 MHz used. It was taken care to use seven nodes per shortest 
wavelength for discretization. It was found that the solutions are 
stable for the resolution using the Newmark’s integration scheme. The 
size of the domain modelled is 30mm x 80 mm in all cases. The details 
of various simulations are given in Table 4.3. 


TABLE 4.3 

FE simulation details for various defects in composites 


Problem 

No. of 

elem. 

No . of 

nodes 

At 

/.j sec 

flaw dimension 

mm 

1 . Drilled Hole 





a. With Hole 

6162 

6332 

0.04 

3.2 

b. Defect Free 

6072 

6230 

0.04 

1 


2. Delamination 




r 

a. With Delamina 

7500 

7691 

0.04 

to 

o 

b. Defect free 

7500 

7676 

0.04 



It was found that for the faster QL wave, the displacements along 
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the fibers are prominent and the wave profile of the same is stretched 
along the fiber direction. For the slower wave QT the displacement 
perpendicular to the fibers are found to be prominent [38] . This is 
the case where the transducer is modelled as a line source with a 
sinusoidal excitation for one full cycle. In the present analysis the 
transducer modelling is as described in eqn. 3.14 and the band width 
of the wave is just equal to the major dimension of the defect 
considered. “ 

4.2.1 Scattering by a hole in a layered composite material 

The wave interaction by a drilled hole in a cross-ply laminate 
made of Graphite/Epoxy has been studied. In the first ply the L-T 
(Longitudinal and Transverse) properties and in the second ply T-T 
properties of the material are used. In the first ply the X-axis is 
along the fiber direction. The hole position is modelled at the 
interface of the layers. 

In the first layer the material is highly anisotropic and the QL 
and QT wave velocities are different in different directions. [Table 
4.2]. The second layer is isotropic in the plane considered and the 
longitudinal and transverse wave velocities are calculated to be 2560 
m/sec and 2130.0 m/sec respectively. 

Figs 4.43 and 4.44 show the X and Y displacement fields when the 
transducer pulse polarised along the propagation direction interacts 
with the hole. The different shapes of the wave profiles can be 
observed in the tw’o layers even at the 6.0 ^sec shot when the wave has 
crossed the hole. A careful observation reveals a sudden change in 
the displacements at the interface. The wave propagation is at an 
angle to the fiber direction in the first ply L-T. In the 8.4 jjsec 
shot, the scattered wave has reached the front wall in the L-T 
laminate where as there is an interference of scattered waves with the 
back wall echos in T-T ply. 

The scattered X and Y displacements were shown in Figs 4.45 and 
4.46. The scattered displacements are highly direction oriented in the 
first ply and are radially distributed in the second layer. Multiple 
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reflections can be observed from the hole surface to the front and 
back walls of the domain. 

Fig 4.47 shows the polar plot of the scattered amplitudes. In 
order to eliminate the scattering from the back and front walls, the 
integrated values oi scattered displacements up to 8 lusec areonly 
considered. When compared with the isotropic case, it contains many 
lobes. The lobes are very narrow in the back-scatter. The scattered 
amplitude is seen to be maximum in the back-scatter i.e. along 180^ 
and the next maximum is likely to occur along 90*^. 

4.2.2 Scattering of SV pulse by a hole in L-T & T-T laminate 

Figs 4.48 and 4.49 show the X and Y displacements when the 
applied pulse is in perpendicular to the propagation direction. The 
X-displacement is predominant in the longitudinal direction in L-T 
lamina and it is seen to be radial in the T-T ply. Even at the 5.0 
^sec shot also one can observe the scattering from the hole. The 
sudden change in the wave fronts is very clear at the interface. The 
Y- displacement profile is very interesting. As such the X 
displacement of the QL wave is seen to be propagating along the 
Y-axis in the L-T lamina. A clear distinction of QL and QT waves is 
not there in these snapshots. 

The scattered X and Y displacements are shown in Figs 4.50 and 
4.51. The scattering is predominant here right from the 6.0 fisec shot, 
i.e. from the first few micro seconds itself, as the X displacement is 
maximum along the fiber direction. This elongation of wave front in a 
particular direction explains the phenomenon of beam skewing in 
ultrasonic testing of anisotropic materials. Here also the scattered 
waves are radially distributed from the hole surface in the second ply 
T-T. 


A polar plot of scattered displacements is plotted in Fig. 4.52. 
It also has contained many lobes like structures. The scattered energy 
is seen to be maximum along 50' direction. Next maximum is likely to 
occur along the 90® direction. 
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A- sacn display at a point near to the hole and at 45*^ to the line 
oi propagatioii (in tlie first ply) is shown in fig. 4.53. When compared 
to the results ot Jalleel [38J , wherein the force pulse was at a 
point, the effect of distributed input force pulse is noticed from 
these plots. 

4.2.3. Scattering by delamination in a layered Composite Material 

The wave interaction by delamination in a cross-ply laminate 
made of Graphite/Epoxy has been studied. The geometry of the ply is 
shown in the fig. 4.54. For the first ply T-L properties of the 
material and for the second ply T-T properties are used. 

The X and Y displacement fields at different times are shown in 
Fig. 4.55 and Fig. 4.56 respectively. The evolution of QL wave into 
first medium is visible from the first two snapshots. The reflection 
of QL wave from the side walls is visible from the third and fourth 
shot. As the transverse direction of both fibers is coinciding with 
the propagation direction, the effect of interface is not seen in 
these plots. It can be noticed from the third shot that there is a 
strong scattering of QT wave from the delamination crack. Also 
observed is the effect of pulse input pulse band width. 

The scattered X and Y displacement fields are shown in Fig. 4.57 
and Fig. 4.58. The scattered Y displacement is slightly radially 
distributed in the T-T ply, whereas the X displacement is predominant 
perpendicular and along the fiber direction in both the plies. There 
is not much change in the distribution of scattered X and Y 
displacements . 

Differential scatter curve for this interaction is as shown in 
the fig 4.59. It is also having lobes like structure. The scattered 
energy is seen to be negligible or zero along the 90® direction. As 
observed in the scattered displacements the scattered energy is 
dominating in both 0® and 180® directions, which can be sensed by the 
transmitter itself in the pulse-echo method. 

A- scan display with and without interface crack at an inside 
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point is shown in fig. 4.60* For the sake of comparision the X-* and Y- 
displacement values with and without the interface crack are shown. 
The direct and the reflected QL wave is marked in the figure. The 
presence of QT wave is not felt in this A- scan. 

4.2.4. Scattering of SV wave from an interlaminar crack in T-'L & L~T 
Laminates 

The scattering of a pulse, which is polarised in a direction 
perpendicular to the propagation direction, when it interacts with an 
interlaminar crack in a T-L and T-T cross ply laminate are as shov/n in 
figs 4,61 and 4.62. Here also the scattered amplitudes are seen to be 
maximum along the interface. 

Fig 4.63 shows the angular distribution of scattered 
displacement amplitudes. The scattered power appears to be predominant 
in 135^ direction. When compared with the earlier case, a considerable 
amount of energy is reflected in this case along 90^ direction. Only 
a little amount of energy is scattered in 0^ and 180^. 


Y 
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30 mm 
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FIG. 4.1 GEOMETRY OF DEFECT-FREE ISOTROPIC SOLID 


I 
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riME * 1-8 M sec e* 3.6 (J sec 



Defect-free Isotropic Solid 




TIME » 3,6 lu sec 



Fig* 4.5 X- Displacement! Rayleigh Wave Interaction with 
a Surface Crack. 







TIME « 3.6 



4.6 Y-Displacement; Rayleigh Wave Interaction with 
a Surface Crack. 



1.8 ^ sec TIME » 3.6 fu sec 



Scattered X-Displacement Field: Rayleigh 
Interaction with A Surface Crack. 









Rayleigh wove 
transmitter / receiver 



Fig. 4.9 Wave Field Components which result from the interaction 
of pulsed Rayleigh Waves with A Surface Crack. 
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TIME *** 3.6 



X-Displacement Field: Compression Wave Interaction 
with A Surface Crack. 






TIME » 1.8 u sec TIME * 3,6 lu 



with A Sui'face Crack. 







time * 3.6 M &ec 



TIME » ^ sr©c TIME * 3.6 ^ 





Scattered Y-Displacement Field: Compression Wave 
Interaction with A Surface Crack. 





Fig. 4.15 (a) Interpretation of Photoelastic Visualisation of the 

Field After A Plane Compression Wave has 
Scattered from a Slit, (h) Contours of Compression 
Wave to correspond with (a).(c) Contours of Shear 
Wave to correspond with (a). 
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Fig. 4.17 Polar Plot of Scattered Amplitudes by A Compression 
Wave Interaction with A Surface Crack. 
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Fig. 4.18 Finite Element Discretization of a Plate with Hole. | 
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TIME m l.e u 



Fig. 4.19 X-Displacement Field: Compression Wave Interaction 
with Hole in Isotropic solid. 



TIME “ 25.B 



Fig. 4.20 Y-Displacement Field: Compression 

with Hole in Isotropic solid. 


1.8 u sec TIME » 2.B ^ sec 



4.21 Snapshots of the Scattered X~ Disp. Field : Compression 
Wave Interaction with Hole in Isotropic Solid* 






TIME * 2.8 fu sec 



Fig. 4.22 Snapshots of the Scattered Y- Disp. Field : Compression 
Wave Interaction with Hole in Isotropic Solid. 








Fig* 4*24 A-SCAN plot at a point before the hole surface 
(Compression wave interaction) 












2.0 mic.sec riME =3.8 mic.sec 



Y-Displacenient Field: Compression Wave Interaction 
with Hole in Isotropic solid. (Resonance Condition) 







.0 mic.sec time =3.8 mic.sec 



Fig. 4.27 Snapshots of the Scattered X- Disp. Field ; Compression 
Wave Interaction with Hole in Isotropic Solid. 





joic.sec TIME — 3.8 mic.se 



Fig. 4.28 Snapshots of the Scattered Y- Disp. Field : Compi-ession 
Wave Interaction with Hole in Isotropic Solid. 

(Resonance Condition) 









1,8 u sec TIME « 2.8 ^ sec 



Figi 4.30 Snapshots of the Scattered X- Disp. Field : Plane Shear 
Wave Interaction with Hole in Isotropic Solid. 
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Fig. 4.31 Snapshots of the Scattered 1 Disp. lield 

Wave Interaction with Hole in Isotropic Solid. 
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Fig. 4.32 Polar Plot of Scattered Amplitudes by A Shear 
wave Intei-action with A Hole. 
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Fig. 4.33 


f A Wedge Crack 

Finite Element Descretization 

(Body Defect). 
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Fig. 4.38 Angular Distribution of Scattered Amplitudes after a 
Plane Compression Wave interaction with Wedge Crack. 











Y-DISPLFKZMENT 
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Fig. 4.42 Angular Distribution of Scattered Amplitudes after 
Plane Shear Wave interaction with Wedge Crack. 






TIME - 6.0 u sec TIME « 7.2 ^ sec 



Fig. 4.43 X-displacement field : Plane Compression wave interaction 
with a hole in Layered Structui-e. 









Scattered X- disp. field : plane Compression wave 
interaction with a hole in Layered Structure. 







7.2 H sec 






TIME - 6.0 (J sec TIME 7.2 fj sec 



Fig. 4,49 Y-displacement field : Plane Shear wave intei’action 





time » 6.0 fJ sec time * 7.2 sec 



Fig. 4.50 Scattered X- disp. field : plane Shear wave 
interaction with a hole in Layered Structure. 


TIME • 6.0 ju see TIME ** 7.2 sec 



Fig, 4.51 Scattered Y- disp. field : plane Shear wave 









iretization of Layered Composite 
Interface Crack. 
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TIME = 2,0 (u sec TIME = 4-0 (LJ sec: 



Fig. 4.55 X- displacement field : Plane Compression wave interaction 
with an Interface Crack in Layered Structure. 









TIME = 2.0 sec TIME = 4.0 f.i sec 



4.57 Scattered X- disp. field : plane Compression wave 
interaction with an Interface Crack in Layered Structure. 


time * ^ TIME = 4.0 AJ sec 



4.58 Scattered Y- disp. field : plane Compression wave 
interaction with an Interface Crack in Layered Structure. 
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A-SCAN PLOT A-SGAN PLOT 




Fig. 4.60 A- Scan plots comparision in a T-L & L-T composite 
laminate with and without interface crack. 



4.0 M sec time = 6.0 /u sec 



4ibl ScattGi'ed X~ disp. field i plane Shear wave 





IME » 4.0 M time = 6.0 M 





4.62 Scattered Y- disp. field : plane Shear wave 
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5,1 Conclusions 

A two-dimensional finite element code has been used to solve the 
hyperbolic equations governing ultrasonic wave motion and to study the 
scattering in isotropic (Aluminium) and orthotropic (Graphite/Epoxy) 
media* A number of simulations have been performed to demonstrate the 
power of the model in describing the scattering phenomenon by surface 
and body defects. The ’’snapshots’’ of the wave field give insight into 
the physical process involved. The model can be used to optimize the 
parameters of a ultrasonic NDT system. 

1 . A clear picture of the scattered waves by flaws is achieved by 
subtracting the incident field from the total field. 

2. The relation of finite element mesh design dispersion phenomenon of 
wave propagation is established and its effect was minimised by 
carefully modelling the flaw free plate domain. 

The following conclusions are drawn from the results of isotropic 
solids : 

1 . The finite element code is capable of predicting the complex mode 
conversion of incident waves by surface cracks and body defects. The 
mode-conversion of Reyleigh wave into a shear and compression waves and 
their propagation behavior have been, exemplified. This understanding of 
mode- conversion of compression wave by a surface crack helps in the 
characterization of surface defect. 


2. The Rayleigh wave is not clearly seen in the scattered waves by the 
body defects, this may be due to the less availablility of surface area. 
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3. An estimate of the total energy scattered by the flaws in different 
directions was made by the differential scattering cross-section. The 
disti’ibut ion is seen to be smooth in case of surface cracks and wedge 
crack. For the case of a hole sharp peak amplitudes were observ'ed in 
both forward and back-scatter fields, when the wavelength is much larger 
than the size of the hole. Resonance effect is observed when the 
wavelength of the input pulse is equal (appro. ) to the hole diameter. 
The effect of wave polarisation on the differential cross-section plots 
was observed. 

A marked difference of the wave scattering in composite material from 
that of an isotropic material is observed, due to highly anisotropic 
nature of composite material. The following observations are made for 
the ultrasonic wave propagation and scattering in orthotropic composite 
materials : 

1. The scattering wave front propagation is clearly seen in the 
different layers { T-T and L-T; T-T and T-L ) . It is highly directional 
dependent. The effect of wave polarisation is also observed. 

2. The surface waves in the orthotropic graphite /epoxy is found to be 
not as prominent as in isotropic media. 

3. The complex lobe structure in the different scatter plots helps in 
detecting the flaw by properly locating the ultrasonic receivers. 

5.2 SCOPE FOR FUTURE WORK 

Future work seems to have a lot of potential, in the following 
directions ; 

3 . A careful study can be done to find out the transducer parameters 
like beam steering angle, beam width variation etc. 

2. In order to optimise the selection parameters like frequency of pulse 
etc a study of the variation of scattering cross sections in a 
particular direction with the variation in the parameter ka can be done. 
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3. Efficient absorbing/nonreflecting boundaries should be found to 
eliminate the spurious reflections from the boundaries. In the 
investigation of these boundaries one should keep in mind the 
computational cost and their ease of incorporation into the model. 

Other wise, it may be advantageous to combine the boundary element 
method with the finite element method in an effort to create an optimum 
hybrid scheme which can be characterised by the FEM solution in the near 
field and BEM solution for the far field. 

4. In the present analysis, the flaws are assumed to be present through 
out the thickness (because of 2-D plane strain approximation), A three 
dimensional Finite Element model apart from allowing to model a variety 
of defects also takes into account their thickness. Italso considers 
the effect of antiplane displacement component w. 

5. The characterisation of the defect/flaw must be done quantitatively 
to assess its severity as it is an ultimate aim of this work. 

6 . Laboratory experiments can be carried out using the Ultrasonic 
scaning equipment in order to validate the results. 
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scaning equipment in order to validate the results. 


otrppDTrMriro 


1. Cantwell, C.J., and Morton, J., "The significance of damage and 
defects and their detection in Composite Materials", J. Strain 
Analysis, Vol.Z7(l), pp zy-42, lyy'Z. 

2. Christensen, R.M,, "Mecljanics of Composite Materials", John Wiley 
Publications, lyTB. 


3. Ki'autkramer , J., and Krautkramer, H., "Ultrasonic testing of 

materials”, Spi inger-verlag , Fourth edition, iy83. 


4. Posakony, G.J., "Engineering aspects of ultrasonic piezoelectric 
transducer design" in IEEE Ultrasonic symposium proceedings, New 
York, pp 1-9, 1975. 


5. Pao , Y'.H., "Elastic waves in solids", J. Appl . Mech.,Tr. 
ASME,Vol.50,pp 1152-1164,1983. 

6. Gubernatis, J.E., "Elastic Waves scattering methods: 
Assesments and suggest ions" , In Rev. of progress in Quantitative 
NDE, 5A, Thompson, D . 0 ., and Chimenti,D.E. ,eds. , Plenum press, Ni,pp 
21-29,1986. 


7. Bond, L.J., Punjani, M. , and Saffari, N ., "Ultrasonic Wave 
propagation and scattering using explicit EDMs", In Math. Model!, 
in NDT , Blakemore, M., and Georgiou, G.A. , eda. , Clarendon press, 
Oxfprd,pp 81-124,1988. 


8. Coffey, J.M. and Chapman , R . K . , Appl icat ion 

Scattering Theory for smooth flat cracks to 
prediction of ultrasonic defect detection and 
Energy, Voi.22, pp 319-331, 1983. 


of Elastic W'ave 
the quantitative 
sizing", Nuclear 


9. Pao, Y.H., and Mow, C.C. , "Diffraction of Elastic and 

dynamic stress concentrations", Crane and Russak press, N\, 19/ 


Gubernatis, J.E., Domany , E 
"The Born approximation in 

Apll 


10 . 

H . , 

elastic waves by flaws , 


, Krumhansl , vi 
the theory of 
Phys. ,Vol.48.pp 


.A., and Huberman, 
the scattering of 
2812-2819,1977. 


11 . 

1 ", 


Pao , Y . H . 
J . Accoust 


and Saches, 
Soc . Am. , VO 


W., "Theory of Elastic Wave 
1.56l5).pp 1478-1486, 1974. 


12. Varatharajulu, V. and 
elastic waves. I Theory , 
556-566, 1976. 


Pao, Y.H., "Scattering 
J. Accoust. Soc. Am., 


Diffraction 

matrix for 
Vol.60, pp 


126 



au. , rv.n., A new way to calculate scattering of acoustic 

and elastic waves 1. Thery illustrated for scalar waves", J. Appl . 

Phy., Vol. 51 ,pp 825 - 834 , 1980 . 

14. Achenbach, J-D., Gautesen, A.K., and Mcmaken, H,"Ray methods 
for waves in elastic solids-with applications to scattering by 
cracks" , Pitman press, 1982 

15, Keller, «1 • B "Diffraction by an apperture " , J , Appl. 

Phys,Vol.28,pp 426-444,1957. 


16, Oglivy, 
austenitic 
Ultrasonics , 


J.A., " Ultrasonic beam profile and propagation in an 
weld using a theoritical ray tracing model", 
Vol.24, pp 337-347, 1985. 


17. Bond, L.J., " Numerical techniques and their use to study 

wave propagation and scattering - A review", in elastic waves and 
ultrasonic NDE, Datta, S.K., Achenbach, J.D., and Rajapakse, Y.s., 
North-Holland , pp 17-27, 1990. 


18. Chin, R.C.Y., Hedstrom, G. , and Thigpen, L. , "Numerical 
Methods in Seismology", J. Comput. Phys . , Vol.54, pp 18-56, 1984. 

19. Manolis, G.D., and Bekos, D.E., " Boundary element methods in 

Elastodynaroics" , Unwin HYman press, London, 1988. 


20. Schafbuch, P.J., Rizzo, F.J., and Thompson, R.B. , " BEM 

solutions for elastic waves scattering in 3D", Int . J. Num. Meth. 
Eng., Vol.36, pp 436-455, 1993. 


21. Hall, W.S., and Robertson, W.H., "BEA of Acoustic wave 
scattering", in math. meth. in NDT, Blakemore, M., Georgious, 
G.A., eds., Clarendon press, pp 357-361, 1988. 


22. 

Alterman , Z . , 

and 

KAral , 

F.C. , 

” Propagation 

of elastic waves 

in 

367 

layered media 
-398, 1968. 

by 

FDM" , 

Bull. 

Seism. 

Soc . 

Ame . , Vol . 58 , 

PP 

23 . 

Boore, D.M. 

II 

> 

FDMs 

for 

seismic 

wave 

propagation 

in 


Heterogeneous materials", in methods in computational physics. 
Bolt, B.A. , ed. , Academic press, NY, 1972. 

24. Blake, R.J., Bond, L.J. and Downie, A.L., "Advances in Numerical 
Studies of Wave Propagation and scattering", in Review of progress in 
QNDE. vol.l., Thompson, D.O. and Dale, E.C., ed., pp 157-164, 1981. 

25. Sochacki, J., et al , "Absorbing Boundary Conditions and 
Surface Waves", Geophysics, Vol.52, pp 60-71, 1987. 

26. Bond, L.., and Saffari, N., "Mode-Conversion Ultrasonic 
Testing", in Research Techniques in NDE, Vol. 7, Ed. R.S. Sharpe, 
pp 145-189, 1984. 


127 



27. Temple, J.A.G., "Modelling the propagation and scattering of 
Elastic waves in inhomogenious and anisotropic media", J. Phys. 
D:Appl. Phys., Vol.21, pp 859-873, 1988. 

28. White, R.M., "Elastic wave scattering at a cylindrical 
discontinuity in a solid", J. Acoust. Soc. Am., Vol.30, p 771-785, 
1958. 


29. Kasaburo Harumi , " Computer simulation of ultrasonics in a 

solid". Materials Evaluation, Vol.44, pp 1086-1110, 1986. 


30. Temple, A., and Ogilvy, J.A., " Numerical techniques for wave 
propagation and scattering in inhomogeneous anisotropic materials", 
in elastic waves and ultrasonic NDE", Datta, S.K., Achenbach, J.D., 
Rajapakse, Y.S., eds., North-Holland , 1990. 


31. Smith, W.D., " The application of FEA to body wave propagation 
problems", Geophys . J. R. Astr. Soc., Vol.42, pp 747-768, 1975. 

32. Ludwig, R., and Lord, W., "A Finite Element Formulation for the 
study of ultrasonic NDT systems", IEEE Trans. on Ultrasonics, 
Ferroelect rics and Frequency control, Vol.35, pp 809-820, 1988. 


33. Luidwig, R. , Moore, R. , and Lord, W., " Transducer models for 

the Finite Element Simulation of ultrasonic NDT phenomena , in 
Review of progress in QNDE, Thompson, D.O., and Chimenti, C.E., 
eds., NY, Plenum press, V 0 I. 6 A, pp 649-655, 1987. 


34. You, Z., and Lord, W. , " Finite 

interaction with cracks, in Review of 
D.O., and Chimenti, D.E., (editors), 
109-116, 1989. 


Element study of Elastic wave 
progress in QNDE", Thompson, 
NY, Plenum press, V 0 I. 8 A, pp 


35. Seron, 
propagation 
pp 359—368 , 


F.J., etal., " Finite Element Method for Elastic wave 
*Communications in applied numerical methods, V 0 I. 6 , 
1990. 


36. You, Z., 
propagation in 
QNDE", Thompson 
1988. 


elal., " Numerical modelling of elastic wave 
anisotropic materials", in Review of progress in 

. D., and Chimenti, C.E., eds., Vol.7A, pp 23-30, 


37. You, Z., and Lord, 
diffraction tomography," 
D.O. , and Chimenti, C.E., 


" A Finite Element test bed for 
in Review of progress in QNDE, Thompson, 
eds., NY, Plenum press, Vol.9, 1991. 


38. Jaleel, K.M.A. 
Element , Modelling 
Materials", M.Tech. 


of Wave Interactions 
Thesis Report, March 


with Cracks in 
1992, IIT Kanpur. 


"Finite 

Composite 


39 Barratt, P.J., and Collins, W.D., "The scattering 
IV. Tooele in’an elastic solid for Plane harmonic 
Camb. Phil. Soc., V 0 I. 6 I, PP 969 981, 


cross section 
waves", Proc . 


128 



40 . Krenk , s . , 
circuiai- crack" 
198;d . 


and Schmidt, H. , "Elastic wave scattering by a 
Phil. Trans. R. Soc . bond., Vol.A308, pp 167-i9B, 


41. Su , J.H., Vardan, V.V., and Varadan, \'.K., "The Unimoment method 
lor elastic wave scattering problems", in Review of progress in 
QNDEVl, Thompson, D.O., and Chimenti, C.E., eds., pp 137-143, 1981. 

48. Goetschel, D.B., Dong, b.b., and Muki , R., "A global local 
linite Element Analysis of axisymmetric scattering of elastic- 
waves", J. Appl . Mech. , Trns. ASME., Vol.49, pp 816-820, 1982.43. 

43-Eari , f. Graff., " Wave motion in elastic solids". Clarendon 
pres s ( Oxford ) , 1975. 

44. Achenbach, J.D., " Wave propagation in solids", North-holland , 

1970 . 

45. Musgrave, M.J.P., "Elastic Waves in Anisotropic Media", in 
Progress in Solid Mechanics Vol . 2, Ed. Sneddon, I.N., and Hill, R. 

46. Mai, A.k., and Singh, S.J., "Deformation of Elastic Solids", 
Prentice Hall Publications, New Jersey, 1991. 

47. Hearmon, R.F.S., "An introduction to Applied Anisotropic- 
Elasticity", Oxford University Press, London, 1961. 

48 Gubernatis, J.E., Domany, E. and Krumhansl , J.A., " Formal 
aspects of theory of the scattering of ultrasound by flaws in 
elastic materials", J. Applied Physics, vol. 48, pp 2804 - 2811, 

1977 . 


49. Newmark, N.M., "A method of computation for structural 
Dynamics", J. Eng. Mech.Div., ASCE, Vol. 85, pp 67-94, 1959. 

50. Smith, W.D., "A nonreflecting plane boundary for wave 
propagation problems", J. comp, phys., Vol. 15, pp 492-503, 1974. 

51. Reynolds, A.C., "Boundary conditions for the numerical solution 
of wave propagation problems". Geophysics, Vol. 43, pp 1099-1110, 
1978 . 

52. Langenberg, K.J. at al , Numerical Modelling of Ultrasonic 
Scattering", in Mathematical Modelling in Non-Destructive Testing, 
Eds. Blakemore, M. and Georgiou , G . A. , Clarendon Press, pp 125-174, 
1988 . 

53. Harker, A.H., "Numerical Modelling of the Scattering of Elastic 
Waves in Plates", J. Nondestructive Evaluation, Vol. 4, No. 2, pp 
89-105, 1984. 


129 



